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We present a theoretical analysis of three-dimensional (3D) matter-wave solitons and their sta- 
bility properties in coupled atomic and molecular Bose-Einstein condensates (BEC). The soliton 
solutions to the mean-field equations are obtained in an approximate analytical form by means 
of a variational approach. We investigate soliton stability within the parameter space described 
by the atom-molecule conversion coupling, atom-atom s-wave scattering, and the bare formation 
energy of the molecular species. In terms of ordinary optics, this is analogous to the process of 
sub/second-harmonic generation in a quadratic non-linear medium modified by a cubic nonlinear- 
ity, together with a phase mismatch term between the fields. While the possibility of formation of 
multidimensional spatio-temporal solitons in pure quadratic media has been theoretically demon- 
strated previously, here we extend this prediction to matter-wave interactions in BEC systems where 
higher-order non-linear processes due to interparticle collisions are unavoidable and may not be ne- 
glected. The stability of the solitons predicted for repulsive atom-atom interactions is investigated 
by direct numerical simulations of the equations of motion in a full 3D lattice. Our analysis also 
leads to a possible technique for demonstrating the ground state of the Schrodinger-Newton and 
related equations that describe Bose-Einstein condensates with non-local inter-particle forces. 

PACS numbers: 03.75.-b, 03.75.Lm, 42.65.Tg 



I. INTRODUCTION 

Recent developments in Bose-Einstein condensation of 
alkali gases include the possibility of coherent molecule 
formation 0,@]j or superchemistry Q via Bose-enhanced 
chemical reaction at ultra-cold temperatures. The rele- 
vant parametric quantum field theories flSSIEB and 
their classical non-linear optical analogs [a, B 13 are 
the subject of much current attention, due to the possi- 
bility of stable, bright, higher-dimensional solitons (also 
referred to as solitary waves) 0, 0, El El • Parametric 
solitons in nonlinear optics with quadratic nonlinearity 
have now been observed experimentally in two transverse 
dimensions, both as spatial and temporal solitons |T^|. 

The dynamical equations for parametric solitons are an 
example of a classically non-integrable field theory which 
generally needs to be treated numerically. The formation 
of three-dimensional (3D) localized solitons is a subject 
of much intrinsic interest in mathematical and non- linear 
physics, and it is intriguing that no integrable models 
supporting them appear to exist. Despite the absence of 
integrability, the quadratically coupled equations of para- 
metric nonlinear optics and coherently coupled atomic- 
molecular Bose-Einstein condensate (BEC) systems ap- 
pear to be the simplest physically relevant Hamiltonian 
models having 3D localized solitons |lj| . 

They also provide an experimental route towards 
demonstration of the closely related ground state of the 
Schrodinger-Newton (SN) equation [14( introduced to de- 
scribe gravitationally bound Bose gases, and later revived 
by Penrose and others [TH as a possible model of 
the collapse of the quantum-mechanical wave-function. 
We also show that there are parallels with more general 
mean-field models of Bose gases having a combination 



of short-distance repulsion and finite-range Yukawa-like 
attractive interactions. These more general models may 
also have astrophysical or quantum-mechanical signifi- 
cance. 

The surprisingly close parallels to nonlinear optics, 
in the related fields of quantum many-body theory and 
atom optics, have now resulted in the emergence of a new 
field of research - nonlinear atom optics with paramet- 
ric nonlinearity. The first step towards seeing molecular 
condensation was recently undertaken in transient exper- 
iments with a Bose-Einstein condensate of 85 Rb atoms 
|17| , in which interference measurements were indicative 
of coherent molecule formation. More recent experiments 
with 133 Cs, 87 Rb, and 23 Na [3, as well as with degen- 
erate Fermi gases of 40 K and 6 Li atoms have pro- 
duced even larger fraction of ultra-cold molecules, as well 
as Bose-Einstein condensates of bosonic molecular dimers 
composed of fermionic atoms. All these experiments have 
employed magnetic Feshbach resonances, which appear 
to be more successful at present than the alternative Ra- 
man photo-association scheme [2(| • 

In addition to this remarkable experimental progress, 
the original effective quantum field theory 0, 0, 0, 01E El 
for coupled atomic-molecular BECs has also been devel- 
oped; it now incorporates renormalization, the treatment 
of intrinsic pair correlations, quantum fluctuations, and 
thermal effects (see, e.g., [pllllllHIllllllHlll 
l29j . as well as a recent review paper [3(J for further ref- 
erences). 

In the near-classical limit of large numbers of atoms or 
photons, the relevant equations for parametric solitons 
are those of mean-field theory, which is a modified ver- 
sion of the Gross-Pitaevskii (GP) equation (in the atomic 
BEC case) or the nonlinear Schrodinger equation (in the 
photonic case) . These are in fact identical equations, ex- 
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cept expressed in the different languages of condensed 
matter physics or photonics. The modification consists 
of the addition of a parametric nonlinear term analogous 
to a quadratic nonlinearity in nonlinear optics. This cou- 
ples the atomic and molecular fields together by means of 
a coherent inter-conversion process. The parametric cou- 
pling acts as an energy- lowering 'glue', which can permit 
stable mutually-trapped BEC solitons to form in 3D. In 
the absence of gravity, this would imply the possibility 
of localized matter waves in free space without external 
trapping potentials. Unlike the usual GP equations, the 
existence of attractive forces in three dimensions does 
not result in a catastrophic collapse, provided the s-wave 
scattering length is non-negative. 

In nonlinear optics, various aspects of the competition 
between quadratic and cubic nonlinearities on soliton for- 
mation have been studied in Refs. |3l| . either in lower 
dimensions (ID or 2D ) or in cases of special relations 
between the system parameters. Matter- wave solitons in 
3D coupled atomic-molecular BECs have been studied in 
Refs. 0, IE 0| and , but only for specific values of 
the s-wave couplings. 

In the present paper, we extend these results to the 
case of arbitrary interaction strengths for the atom- 
molecule coupling and for the repulsive atom-atom s- 
wave scattering, as well as for arbitrary energy mismatch 
between the atoms and molecules. We analyze the su- 
perchemistry equations in the mean- field limit, to obtain 
the precise conditions under which 3D atomic-molecular 
BEC solitons can form. Approximate soliton solutions 
are found analytically, by means of a variational approach 
with a Gaussian and an exponential ansatz. We then nu- 
merically study the dynamical stability of the resulting 
solitons on a 3D lattice, together with comparing the 
results with exact numerically solutions. 

We find that there are large regions of stability in pa- 
rameter space, depending on the energy difference be- 
tween the atomic and molecular condensates, the num- 
bers of atoms involved, and the coupling strengths. For 
simplicity, the analysis only includes repulsive s-wave 
scattering between the atoms, and assumes no other s- 
wave interactions. While more general s -wave interac- 
tions are simple enough to include, we have focused on a 
relatively simple case here in the interest of keeping the 
parameter-space manageable. 



II. THE MODEL 

We start with an effective field theory model for a co- 
herently coupled atomic-molecular system pj . The model 
and the obtained results can easily be adopted to describe 
certain cases of non-linear optical interactions of second- 
and sub- harmonic waves in a nonlinear crystal 0, El • In 
the atom-molecular case, the model refers to a type of su- 
perchemistry in which an atomic condensate is able 
to coherently and reversibly inter-convert with a conden- 
sate of diatomic molecules. 



There are several possible experimental routes for pro- 
viding this type of coherent coupling, including a Fes- 
hbach resonance (employing a tuned external DC mag- 
netic field), Raman photo- association (involving two ex- 
ternal lasers with a well-defined frequency difference), 
and direct single-photon photo-association (this would 
require an external microwave or infra-red field) HUH 
IllHlllllllilillilllllillillilll^. The first two 
cases have been experimentally demonstrated 0, l20| . 
although not yet the last. In practical terms, the main 
limitation of the model and corresponding experiments is 
the need to minimize incoherent processes like inelastic 
collisions and other loss processes. 



A. Hamiltonian 

In the model, we suppose that each condensate has the 
usual kinetic energy term, atom-atom s-wave scattering 
interactions, and a number-conserving coherent coupling 
of the form $^<&^ , I', where $ represents the atomic field, 
and is the field operator for the molecular dimers. In 
D (D = 1, 2, 3) spatial dimensions, this leads to a model 
Hamiltonian of the following form: 



H 



h 2 h 2 

" 'V$(x)| 2 + - |V*(x) 



2mi 



2m 2 



+l/** t (x)*(x) + — ii$ t (x)l> t (x)l'(x)<l(x) 



^ ($ t (x)$ t (x)*(x) + ^ t (x)$(x)$(x; 



(1) 



Here, mi and 7712 are the masses of the atoms and 
molecules, respectively, V* is the internal molecular en- 
ergy relative to free atoms, and the coupling \ (which 
we assume is positive) describes coherent conversion of 
pairs of atoms into diatomic molecules and vice versa. 
The atomic self-interaction strength, K\x, is proportional 
to the s-wave scattering length. For example, in 3D, 
Kn = ATrhan/mi, where an is the atom-atom scatter- 
ing length which is assumed positive, as is usually needed 
to form a stable BEC in the first place. 

To allow comparisons with the Schrodinger-Newton 
(SN) 0, 0, 0] equation, we will also consider a re- 
lated model in which the interaction term — hx[^^^' + 
* t $<I>]/2 is replaced by - fi-x^^I* + * f ]/2. This mod- 
els a Bose-Einstein condensate ^(x) with short-range in- 
teractions, together with a long-range attractive force 
caused by the exchange of a meson-like particle ^(x). 
In the mean-field theory limit, we call this model the 
Gross-Pitaevskii- Yukawa or GPY model, which is more 
general than the Schrodinger-Newton model. 

In astrophysical situations, the GPY mean-field the- 
ory reduces to the SN model in the combined limit of 
zero s-wave scattering, and an infinitely long-range grav- 
itational interaction with m<i — > 0, leading to an inverse- 
square law |14j . The presence of a short-range interaction 
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s-wave scattering term makes the GPY model more re- 
alistic than the usual SN model. The SN model is used 
to describe a degenerate Bose gas with gravitational self- 
interaction, and has also been suggested as a possible 
mechanism for wave-packet collapse in quantum mechan- 
ics In the one-dimensional case, the GPY theory 
is similar to the nonlinear interactions in an optical fiber 
caused by couplings of photons to phonons [33. At the 
same time, the GPY model allows one to investigate at- 
tractive interactions with more general behavior than a 
simple inverse-square law. 

At the quantum field level, either model Hamiltonian 
implicitly involves a delta-function effective interaction 
between the atoms, and so requires the use of a mo- 
mentum cutoff k max <C 1/fln for selfconsistency. More 
rigorous regularization consists of renormalization of the 
theory [2(j as k max — » 00. In the present paper, however, 
we employ a mean-field approximation, and restrict our- 
selves to the study of solitons that have spatial widths 
much larger than k^ ax so that the cutoff dependencies 
are negligible, while the relevant parameters are the ob- 
served couplings. The mean-field theory is also a high- 
density approximation, since quantum fluctuations and 
correlations are expected to cause quite different ground 
state to appear at low density 0, EU EH • 

A more complete model Hamiltonian should also in- 
corporate atom-molecule and molecule-molecule s-wave 
scattering interactions. However, these greatly compli- 
cate the analysis, without adding much qualitatively new 
physics to the 3D soliton properties studied here. In ad- 
dition, the respective scattering lengths are not known 
yet in most cases. For this reason, we assume that the 
atom-atom scattering is the dominant s-wave interaction 
and simply omit all other s-wave scattering processes, in 
the interest of simplicity. 



B. Mean-field equations 

The corresponding equations of motion for the mean 
fields in the atom-molecular model, following from the 
Hamiltonian and valid at high densities, are 



for the related GPY equations. These have the structure: 



■ fl0(x,t) 

'' dt 



2mi 



where hAu) — is the energy mismatch on converting 
atoms to molecules, and m 2 = 2m \. 

In this model, the total number of particles N (i.e. the 
total number of atomic particles, including pairs of atoms 
inside the diatomic molecules) here is conserved: 

N = N 1 +2N 2 = y"d D x[|0(x,t)| 2 + 2|V(x,t)| 2 ] . (3) 

For completeness, we include the mean-field equations 



■ fl#x,t) 

1 dt 

■ frfr(x,t) 



2m 



-V^-|^ + V*)+«ii|0|V, 



(4) 



In the limit of m 2 — > 0, x ~~ * 00 1 so that Awm2 — * 
and X 2?7l 2 — 47rGm 2 , and assuming that ip is real, we 
introduce a gravitational field potential V g — — x^P- ln 
this long-range force limit, one can apply an adiabatic 
approximation to the second equation, which leads to 
the Poisson equation: 



ih 



dt 



2mi 



V 2 x cl> + V g <t> + hKn\4>\ 2 <l>, 



V 2 Vg = AnGm 



2 H 2 - 



(5) 



This can be recognized as the mean-field equation for 
a BEC having an additional self-gravitational force with 
the gravitational potential energy V g and gravitational 
constant G, as well as the usual GP short-range inter- 
action. Here, the conserved particle number is given by 
N = f g? 3 x \(j)(x, t)\ 2 . In the additional limit of K\\ — ► 0, 
the equations correspond to the time-dependent version 
of the SN equation (l4L fl5L fl^ | in which there is no short- 
range self-interaction. 



C. Dimensionless variables 

1. Atom-molecular system 

In general, atom-molecular solitons may exist with pe- 
riodic phase and frequency uj, so that id<j)/dt = uj<f> and 
idifj/dt = 2uip. 

We introduce a characteristi c time sca le to and a char- 
acteristic length scale do = y/hto/2mi. We also trans- 
form to dimensionless time and position variables, 



T = t/t , 

£i = Xi/do, 



(6) 



and dimensionless fields, 



(7) 



This gives the corresponding equations of motion in di- 
mensionless form, with no reduction in parameter space: 



du 

'j- = 
or 



— V^w + j u u — u*v + an \u\ 2 u, 



.dv 1_ 9 7 I 9 
l d-r--2^ V+ 2 V ~2 U 



(8) 
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Here we have introduced new dimensionless parame- 
ters according to 



a n 

In 

7 



x 2 V 

-wt , 
(2Aw 



(9) 



To make the scaling definite, we can set 7„ = 1 with 
no loss in generality, provided u> < 0. This corresponds 
to a localized bound state with negative energy E = huj; 
we do not investigate the unbound solutions here. The 
choice 7u = 1 also gives a simple relationship between 
the dimensionless parameter 7 and the detuning Au), 



7 = 4 + 2Awt , 



(10) 



corresponding to a shifted energy mismatch. 

For an — (and with an additional scaling of u — > 
uA/2), Eqs. (jHJ are equivalent to Eqs. (1) and (2) of Ref. 
|10| , with the value of the coefficient 8 = 1 . This corre- 
sponds to optical parametric interaction in a quadrat- 
ically nonlinear medium, in which the dispersion coef- 
ficients for the fundamental and second-harmonic fields 
are equal to each other. 

We seek stationary solutions (du/dr = dv/dr = 0) to 
the equations of motion JHJ), i.e. those that have 

u*v + an\u\ 2 u, 

„,2 



V 2 u 



(11) 



These correspond to extrema of the following dimension- 
less atom-molecular Hamiltonian 



H (u,v) = J d D^ 



IVul 



1, 



IVul' 



7, ,2 
2 H 



„_(( U *)2 U + C . C .) + 



an 1 ,4 
-2>l 



(12) 



where the expression for the original Hamiltonian energy 
in terms of the dimensionless variables is 



H 



ft 



D/2 



(13) 



2mit J x" 

A conserved quantity for the (u, v) system, which is 
proportional to the total number of particles N, is 



N'= / d D Z[\u(Z,t)\ 2 + 2\v(Z,t)f 



(14) 



2. Schrddinger-Newton system 

Similarly, stationary solutions to the time-dependent 
SN equation (J5J may also exist with frequency to, so that 
id(j>/dt = uxj), with V g following adiabatically. This trans- 
lates Eq. JSJ directly to the time-independent SN equa- 
tion appearing in Refs. 0,0], for (j) s = (j>exp(iujt): 

ft 2 _o 



2mi 



V 2 K 



AnGm\\<j> s 



(15) 



where E — fku. 

By introducing characteristic time and length scales, 
as in Eqs. ©, and transforming to dimensionless fields 



u = to \/2-KGrn\<\) s , 



ft 9 ' 



(16) 



we obtain the dimensionless time-independent SN equa- 
tions: 



Vjrli = U — VU, 

Vlv = -u 2 . 



(17) 



The normalization is now given by J d D ^u 2 — |£| -1 / 2 . 
Here £ = E/Eq is the dimensionless energy (£ < 0), and 
the energy scale E is defined via Eq — 327r 2 mf N 2 G 2 /h 2 . 

Note that Eqs. l(T7|) are identical to Eqs. l(TT|) for real 
stationary solutions to the atom-molecular system when 
7 = an = 0. 



III. VARIATIONAL ANALYSIS 
A. Gaussian variational ansatz (GVA) 

To analyze the localized soliton solutions to Eqs. JHJ, 
where j u = 1, we introduce an approximate Gaussian 
variational ansatz (GVA). This permits an analytic treat- 
ment of the problem of minimizing the Hamiltonian en- 
ergy. The stability of the GVA solutions will be checked 
numerically by dynamical evolution of the equations of 
motion where the GVA serves as the initial condition. 

The GVA solutions are introduced according to 



u(£,t - 0) = Ae~ a Z , 
v(£,t = 0) = Be- bi \ 



(18) 



where £ = |£|. Here, the parameters a and b must both 
be real and positive for localized solitons, and we assume 
that both the amplitudes A and B are also real and pos- 
itive. This choice in fact already takes care of the op- 
timum relative phase between the atomic and molecular 
fields, where we can without loss of generality take B 
to be real, while the optimum relative phase will dictate 
the phase of A. This immediately leads us to the conclu- 
sion that A 2 must be real and positive too, in order that 
the atom-molecule interaction term remains negative (for 
positive x as assumed here) and permits a minimum in 
the Hamiltonian energy. The sign of A is in fact irrele- 
vant, since the corresponding equations of motion (jSJ) are 
invariant under the sign change of A. 

Substituting the GVA into Eqs. O and O for N' 
and the Hamiltonian energy H^ u ' v \ and taking the in- 
tegrals we obtain, in D = 1,2 or 3 space dimensions: 



D/2 



A 2 



2B 2 



jD/2 b D/2 



(19) 
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2A 2 



B 2 



2A 2 



_ a D/2-l b D/2-l Da D/2 
2 1 + D/2 A 2 B anA 4 



2 

7B 2 

Db°/ 2 ~ D{2a + b) D / 2 ^ D(2a) D / 2 



(20) 



Minimizing H^ u ' v ^ with respect to a, b, A, and B (for 
jiven an and 7) gives the following solution: 



2(£>o + 7 ) 
(£>- 2)6 + 7 



(21) 



an (2a + b) 1+D (Db + 7) [(4 - D)a - 1] 
- 2 1+£, (2a6) D/2 [4a 2 - (D - 2)ab - b] = , (22) 



B = 



{2a + b) D / 2 {Da + l) 



1 - 



A 2 = 



{2a) D / 2 

a 11 {2a + b) D {Db + 1 ) 
2 D {2ab) D / 2 



{2a + b) D l 2 {Db + 1 ) 
{2b) D / 2 



B, 



(23) 



(24) 



In Eq. Ij22(l . the parameter a is to be substituted using 
Eq. ijinjl. so that Eq. l(2^|) (to be solved first) reads as a 
polynomial equation with respect to 6, for given values 
of a 11 and 7. Alternatively, b can be regarded as a free 
parameter and Eq. 11221) be viewed and easily solved with 
respect to an, for given 7 and b. In doing so, only pos- 
itive 6-values have to be considered for physically mean- 
ingful (localized) soliton solutions. 

For certain values of an and 7 this system has a unique 
solution (see Sec. IV A), giving the soliton parameters 
A, B, a and b. The soliton parameters give in turn the 
resulting value of the conserved quantity A/"', Eq. (JTJJJ, 
and the Hamiltonian energy 1)21) \ . 



B. Exponential variational ansatz (EVA) 

It can be shown that any localized stationary solution 
to the equations of motion JSJ (with j u = 1) must possess 
tails decaying according to: 



ii(C»l,r)«e-Ve, 
W (£»l,T)oce-v^/£, 



(25) 



where £ = |£|. This result can be obtained by neglecting 
all non-linear terms in Eqs. (JBJ) at large £, and solving the 
resulting decoupled equations for the stationary states. 

Due to the singularity at origin, direct employment 
of this variational trial function would be problematic. 
However, since it indicates that the soliton tails should 
decay slower than those of the Gaussian trial functions, 



we are motivated to also consider an alternative exponen- 
tial variational ansatz (EVA): 



u (£,r = 0) =Pe- p v^, 
v(Z,T = 0)=Qe-«V^. 



(26) 



As in the case of the GVA, here too we assume that 
p, q, P and Q are all real and positive. For analytic 
simplicity, we let e be an infinitely small length scale, 
which is formally included to ensure that u and v are 
differentiable at £ = 0. We then proceed and evaluate 
the integrals in Af' and H^ u ' v ^ to find that, as e — > 0: 



m' = k l 
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p D 



2tf 



(27) 



K 



D 



P 2 



Q 2 



p 2 



Q 2 



v 

2 D P 2 Q 



D-2 



2q D - 2 
p4 



p D + 1 2q D 
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D 
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2 D+l p D 



(28) 



Here, Ko = 1, 7r/2, it for D = 1, 2, 3, respectively. 

Variational stationary points are then given by the so- 
lution to the following set of algebraic equations: 



1 



2^(g 2 +7) 

(D - 2)q 2 + D 1 



1 



(29) 



2 3D+x (pq) D { [(£> - 2)p 2 + D] (2p + q) - 2Bp(j? + 1)} 
- a n (2p + q) 2D+1 (q 2 + 7) [D - (4 - D)p 2 } = 0, (30) 



{{D-2)p 2 + D)(2p + q) D+1 



1 — ax; 



D(2p) D + 1 

{q 2 + 1 )(2p + q) 2D +^ - 1 
23D+2pD+lqD 



P 



2 _ (q 2 +1 ){2p + q) D 



Q- 



(31) 



(32) 



To compare the EVA solutions with those of the GVA, 
it is necessary to ensure that both have identical J\f' for a 
given pair of the parameters an and 7. As the solutions 
given above result from (unconstrained) variational mini- 
mization with respect to all parameters, this requirement 
will not in general be met. We thus perform a constrained 
minimization with respect to p, q and Q, leaving P to be 
fixed by M' of the associated GVA solution |34|. The 
constrained EVA solution for D = 3 is then given by: 
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-anpF 



48 W Q 

(2 P + q y 
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(33) 
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FIG. 1: GVA solutions existence domain in the (7, an) plane, 
for 7 > 0. The dotted line gives the upper bound on an 
corresponding to the boundary an < I/7. The set of vertical 
lines in the region 7 < 4 (7 = 0.5 - full line, 7 = 2- dashed 
line, and 7 = 3.8 - dash-dotted line), and in the region 7 > 4 
(7 = 4.2 - full line, 7 = 6- dashed line, and 7 = 8- dash- 
dotted line) are to serve for mapping purposes as discussed in 
the text and explained in the captions to subsequent figures. 
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FIG. 2: GVA solution widths a a = 1/V2a and a b = 1/V26 
as a function of an, for different values of 7 corresponding to 
different vertical lines in Fig. Q 
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where we have defined F = Af' /n — 2Q 2 /q 3 . 

We note that a similar exponential variational solu- 
tion for the 'atomic' u -field (though not for the 'molecu- 
lar' v-field) has been used previously for the Schrodinger- 
Newton equation [i"rj| . 



IV. 3D SOLITON PROPERTIES 

A. Existence and properties of GVA solutions 

In order to analyze the properties of the GVA solutions, 
Eqs. H21 |) -l|24 |) . for a given pair of 7 and an, we first note 
that our analysis is restricted to the case of an > 0, 
i.e. repulsive atom-atom interactions including a non- 
interacting limit of an = 0. In addition, we restrict 
ourselves to the three-dimensional case (D = 3) only. 

Next, any localized physical solutions require a and b to 
be both real and positive. The existence of the minimum 
for the Hamiltonian energy, Eq. I|2(J|I . requires that the 
product A 2 B is real and positive too, so that the atom- 
molecule conversion term gives a negative contribution 



to the energy. As we mentioned earlier, this is achieved 
by taking both A and B to be positive. 

To investigate the consequences of these requirements 
in terms of the soliton existence domain in the parameter 
space (an, 7), one has to start from solving numerically 
the polynomial Eq. (|22|) . However, a simpler route that 
allows us to obtain analytical results is to view 6 as a free 
positive-valued parameter and solve Eq. I)22|) for an in 
terms of b and 7. In this case, the GVA solutions can be 
rewritten (for D = 3) in a simpler form: 
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b(5fr + 7) 
2(6 + 7) ' 
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(36) 
(37) 
(38) 
(39) 



where we have substituted the solution for an into the 
expression for B, and the parameter a in Eqs. H37fH39|) 
is to be substituted using Eq. 

In Appendix A, we analyze the above set of equa- 
tions for possible solutions. Restricting ourselves to the 
physically interesting subspace of 7 > 0, we find that 
the soliton existence domain for the GVA is given by 
< an < 1/7- The parameter space identifying this in 
the (7, an)-plane is shown in Fig. ^ Here, the vertical 
lines at 7 = 0.5, 7 = 2, 7 = 3.8, 7 = 4.2, 7 = 6, and 
7 = 8 serve as test lines for mapping purposes discussed 
in subsequent sections. The dotted line gives the upper 
bound on an corresponding to the boundary an < I/7. 
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FIG. 3: Fraction of the number of particles in the atomic field 
Ma/M' for GVA solutions as a function of an, along the lines 
of fixed 7- values as shown in Fig. Q 




FIG. 4: Fraction of the number of particles in the atomic field 
N'a/N' for the EVA solutions as a function of an, along the 
lines of fixed 7- values as shown in Fig. Q . 



the following form (for D = 3) 



Figure [5] represents the GVA soliton widths a a = 
\j\[2a and o\, = 1/V26 as a function of an (0 < an < 
I/7), for different values of 7. Similarly, Fig. ^represents 
the fraction of the number of particles in the atomic field 
M' a /N' versus an, where M' a = J d 3 ^u 2 is proportional 
to the total number of atoms, while Af' = J d 3 £ (u 2 + 2v 2 ) 
is proportional to the total number of atomic particles in- 
cluding pairs of atoms in the molecular component. Due 
to the conserved total particle number, the fraction of 
molecules is found from M' m JM' = 0.5(1 - M'JM'). 

As we can see, for large negative detuning Aco, cor- 
responding to 7 <C 4 (7 > 0), and for vanishing atom- 
atom repulsion (an ~ 0), the atomic fraction is relatively 
small and increases monotonically with increasing 7. In 
all cases, the atomic fraction decreases rapidly as an 
increases, due to the increased energy penalty resulting 
from interatomic interactions. The graphs for the soli- 
ton widths show that the atomic density profiles are in 
general wider than the corresponding molecular density 
profiles, and that the atomic component becomes wider 
and lower in the amplitude in the limit of strong inter- 
atomic repulsion, an — > 1/7. In this limit, it is energeti- 
cally preferable for atom pairs to populate the molecular 
component so that the stable configuration of the sys- 
tem is a pure molecular condensate. We note that this 
is a consequence of the fact that our model neglects the 
molecule-molecule self-interaction completely. 



P = o 



5g 2 + 37 
3 7 



an 



q 

2 V q* 
2 w p 3 q 3 Up 3 



q(p 2 + 3)] 



Q = 



p 2 



(2p + qy(q 2 
(2p + q)\p 2 



7)(p 2 
3) 



3)- 



3-2V(<Z-2p) ' 
{2p + q) 3 {q 2 + 1 ) , 



2 3 q 3 



(40) 
(41) 
(42) 
(43) 



Here, the expression Q41JI for an is obtained from Eq. 
(|30|l and has been further substituted into Eq. I|31l) to 
obtain Eq. fl2*)). The parameter p in Eqs. l|lT |l -l|l3* ]l is 

to be substituted using Eq. I|40p. By treating q as a free 
parameter instead of an (assuming q > for localized 
solutions), we can first solve for p in terms of two inde- 
pendent parameters q and 7, and then proceed to find 
the remaining parameters, an, Q, and P. 

In Appendix B, we analyze the above set of algebraic 
equations and conclude that for 7 > the existence do- 
main for the EVA solutions is identical to that of the 
GVA solutions, and is given by < an < 1/7. 

Figure 01 shows the dependence of the atomic number 
fraction M' a /N' for the EVA solution as a function of 
an, for various 7. As we see, the salient features of these 
curves, as well as the behavior of the EVA soliton widths, 
are similar to those of the GVA solutions discussed in the 
previous subsection. 



V. DYNAMICAL STABILITY 



B. Existence and properties of EVA solutions 



Almost identical arguments relating to the soliton ex- 
istence domain can be made for the EVA solutions. We 
again simplify the analysis by rewriting Eqs. I|29|) - (|32|) in 



In order for these variational solutions to prove useful, 
it is necessary to identify some correlation between their 
dynamical behavior and the existence of actual (exact) 
stable soliton solutions. To this end, we have identified 
the exact stationary solutions numerically, by means of 
the numerical relaxation method, and have checked their 
stability under dynamical evolution for a large number 
of test points in the (011,-7) parameter space within the 
GVA/EVA solutions existence domain, < an < 1/7 
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(a) 




FIG. 5: Example of stable dynamical evolution of the GVA 
soliton. Shown are the particle number densities for the 
atomic (a) and molecular (b) fields for an =0.1 and 7=1, 
with the GVA solution taken as the initial condition. 

for 7 > 0. What follows is an account of the dynamical 
behavior of the variational approximations, together with 
comparisons between this behavior and the existence of 
stable stationary points. 

A. Stability of the GVA solutions 

We have conducted a numerical analysis of the dynam- 
ics of the GVA solution for various (an, 7) pairs lying 
within the existence domain < an < 1/7, for 7 > 0. 
The details of this analysis are as follows. 

We first used Eqs. (|36fl - 139fl to obtain the parameters 
characterizing the GVA solution for each (an, 7) pair in 
question. These solutions were then used, in conjunction 
with the dimensionless equations of motion (jSJ), to form 
a set of initial value problems. The dynamical behavior 
of each Gaussian solution was then determined through 
numerical integration using a spherically symmetric |35j 
semi-implicit algorithm. 

In Figures0through[7|we show typical examples of the 
dynamical evolution of the atomic and molecular fields. 

Figure [S] summarizes the results of our dynamical sta- 
bility analysis, applied to many (an, 7) pairs satisfying 
an > 0.01, 7 > 0.01 and an < I/7. Here, the points 



(a) 




FIG. 6: Dynamical evolution of the atomic (a) and molecular 
(b) field densities for an = 0.01 and 7 = 0.01. This is an ex- 
ample representing strongly 'oscillatory' behavior of the GVA 
solution. 



0.8 T 




FIG. 7: Dynamical evolution of the GVA solution for an = 3 
and 7 = 0.01 representing an example of unstable behavior. 
Shown is the particle number density for the atomic field, with 
similar behavior observed for the molecular field. 



marked with squares, circles or stars represent dynamics 
of the GVA solutions, which have been classified as sta- 
ble, 'oscillatory', or unstable in nature. (This necessarily 
involves a certain degree of ambiguity when distinguish- 
ing between the stable and oscillatory cases.) Here, the 
term oscillatory is used in a broad sense, and does not 
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FIG. 9: Exact soliton dynamics with the exact solitary solu- 
tion (found numerically) as the initial condition, for an = 0.1, 
7=1, and M 1 being fixed to the same value as in the respec- 
tive GVA solution. Shown is the particle number density for 
the atomic field, with similar behavior found for the molecular 
field. 



B. Existence of numerical exact solutions 



10" 2 10° 10 2 

Y 

FIG. 8: Summary of the dynamical behavior of the GVA so- 
lutions for different (011,7) pairs. The squares, circles and 
crosses indicate stable, strongly oscillatory and unstable be- 
havior, respectively. For discussion of the shaded region with 
Eeva < Egva, see text in Sec. V C. 



mean to imply true periodic oscillations around the orig- 
inal GVA solution or the exact stationary solution. The 
term unstable, on the other hand, refers to derealization 
of the GVA solution over short time-scales j^fj • 

Remarkably, the GVA solutions display primarily sta- 
ble dynamics within the < an < I/7 parameter space. 
Exceptions to this are two regions close to the 7 = 0.01 
axis. The lower of these regions (the shaded region in Fig. 
[H] containing circles) contains GVA solutions which, al- 
though remaining localized, display highly oscillatory dy- 
namics. Here, the EVA solutions have lower energy and 
give a better approximation to the exact solutions (see 
Sec. V C for further discussion). The other region con- 
tains unstable GVA solutions (marked by crosses) which 
delocalize rapidly under dynamical evolution. 

We point out, however, that these regions are rather 
small in the physical parameter space (notice the loga- 
rithmic scale in Fig. |8J. The most interesting area in this 
sense is about the 7 = 4 axis, corresponding to a mini- 
mum energy mismatch between the atomic and molecu- 
lar fields, HAuj ~ 0. For example, even for the detunings 
as large as Aw ~ —10 s _1 , which gives an energy mis- 
match comparable in magnitude with typical mean-field 
energies in atomic Bose-Einstein condensates, the corre- 
sponding value of 7 is of the order of 7 ~ 2, for typical 
values of N ~ 10 5 and \ ~ 10~ 6 m 3/2 /s (see also Sec. 
VI). In this physically interesting region, the GVA solu- 
tions are a good approximation to the exact solitons and 
maintain excellent dynamical stability. 



By using the approximate GVA solution as an initial 
guess in the numerical relaxation algorithm, we have in- 
vestigated the shape of the exact stationary solution hav- 
ing the same particle number as the GVA. This constraint 
is used to ensure that the numerically found exact solu- 
tion corresponds to the same set of physical parameters 
as the GVA. The stability of each exact solution was de- 
termined in the same manner as that of the GVA, i.e. via 
real-time dynamical evolution governed by Eqs. JHJ- In 
all cases where a stationary solution was identified but 
found to be unstable, a modified initial guess was found 
that converged to a stable stationary solution. The mod- 
ified Gaussian used in these cases was typically narrower 
and of higher peak density, while having the same total 
particle number. In a small subset of cases, no exact 
stationary solution was obtained (see below). 

FigureEHillustrates the time evolution of one such solu- 
tion. This demonstrates the stability of the exact soliton 
solution, in contrast to those corresponding to energy 
maxima which become delocalized due to the buildup of 
small numerical inaccuracies such as rounding errors. 

In Fig. we summarize the results of the stability 
analysis for the exact solutions and compare them with 
those of the GVA. Stable exact soliton solutions were 
identified for all (an, 7) pairs used in the GVA dynamical 
analysis (see Fig. |SJ, except those lying within the shaded 
region. The boundary of this region was found using an 
extra set of test points for higher accuracy. Distorting 
the initial guess input to the relaxation algorithm did 
not help in finding exact stable solutions in the shaded 
region. Note that, apart from a small number of excep- 
tions near the boundary, all unstable GVA solutions in 
Fig. |H1 which delocalized under dynamical evolution (as 
for the example in Fig.Q) are contained within the shaded 
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FIG. 10: (a) Comparison of stable soliton existence domain 
with the dynamical behavior of GVA solutions. Stable ex- 
act solutions were found for all marked points outside of the 
shaded region. The precise boundary of the shaded region, 
where no exact solutions were found, was identified by ap- 
proaching it from below and from above along the vertical 
lines shown in (b). The lines themselves consist of points, of 
spacing 0.01 in an, for which the exact solutions were iden- 
tified, while the interrupted part of each line corresponds to 
having no exact solutions. 



region in Fig. ^] The physical origin of this instability 
is yet to be understood. 

Thus, one can make the statement that the existence 
and dynamical stability of the GVA solution is strongly 
indicative of the existence of a true stationary soliton 
solution. This applies even to the case of strongly oscil- 
latory GVA dynamics, in the sense that we were able to 
find a stable exact soliton solution whenever the oscil- 
latory behavior of the GVA persisted for long evolution 
time. 



C. Oscillatory dynamics and EVA solutions 

Stationary solitons for (an, 7) pairs in the lower left 
corner of Figure[H]are poorly approximated by the Gaus- 
sian ansatz solutions - a fact revealed by the oscillatory 
GVA dynamics prevalent in this area. Examining the 
profiles of the corresponding numerically-obtained exact 
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FIG. 11: Comparison of the density profiles for the atomic (a) 
and the molecular (b) fields described by the GVA, EVA and 
the exact stationary solutions, for an = 0.01 and 7 = 0.01. 
Note the dramatic failure of the GVA to approximate the 
exact solution in this case. 



solutions suggests that the EVA solutions may provide a 
better approximation in this region. 

In order to test this for a given (an, 7) pair, we need 
to ensure that the GVA and EVA solutions in question 
are being compared for the same value of the parameter 
J\f' corresponding to the total number of particles. Thus, 
we first identified the value of Af' for the GVA solution, 
and then solved Eqs. (|33|) - (|35fl for the parameters of the 
corresponding constrained EVA solution. 

Figure ITU illustrates the improvement in the fit of the 
profile of the EVA solution to that of the exact station- 
ary solution, for an = 7 = 0.01. The corresponding 
GVA solution is also shown for comparison. The result- 
ing reduced amplitude of oscillation in the dynamics of 
the EVA is shown in Fig. [JJJ 

In order to understand and quantify this improve- 
ment, the total Hamiltonian energy corresponding to 
both ansatz solutions has been calculated for a collection 
of points spanning the parameter space under consider- 
ation. The shaded area with Eeva < Eqva in Fig. |S| 
represents the region of parameter space where the con- 
strained EVA solutions have lower energy than that of 
the GVA. This analysis shows that the EVA indeed pro- 
vides a better approximation to the exact solitons in cases 
when the dynamics of the GVA is strongly oscillatory. 

Previous investigations of the SN equation [l^ , which 
has a stationary solution exactly equivalent to ours with 
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(a) 




FIG. 12: Dynamical evolution of the EVA solution for an = 
0.01 and 7 = 0.01. Shown are the particle number densities 
for the atomic (a) and molecular (b) fields. This figure can be 
compared with the strongly oscillatory dynamics of the GVA 
solution in Fig. |S| and shows improvement in the stability of 
the EVA solution due to its lower energy. 



7 = an — 0, have come to the same conclusion that 
over a class of trial functions the exponential ansatz for 
the 'atomic' field provides the best upper bound to the 
ground state energy. The fact that the upper bound pro- 
vided by our solution, — 0.108m 5 G 2 N 2 /H 2 , is higher than 
the value -Q.146m 5 G 2 N 2 /H 2 quoted in is due to the 
fact that we have used variational solutions for both u 
and u-fields, rather than just for u. These upper bounds 
can be contrasted with the exact ground state energy of 
— 0.163m 5 G 2 N 2 /h 2 0] which we have verified using our 
numerical relaxation code. 

The SN helps to also understand the dramatic failure 
of the GVA solutions to approximate the exact solitons 
in the region of small 7 and a. Here, the tails of the 
density distributions become increasingly important, and 
according to Eq. (|^5l) the w-field with 7^0 decays much 
slower than the tail of any Gaussian. In terms of the 
SN equation, this corresponds to evolution of the wave- 
function in a very long-range potential v. 
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FIG. 13: Energy surface topography as one approaches the 
unstable region boundary from below along the 7 = 0.01 line. 
Shown is the dimensionless Hamiltonian energy vs the dis- 
tance s in state space |37|. Different curves correspond (from 
bottom in increasing order) to cm = 0.2, 0.5, 0.7, 1. 



D. Energy surface topography 

In order to better understand the nature of the insta- 
bility marked by the shaded region in Fig. E3 it would be 
instructive to examine the topographical variations in the 
energy surface while approaching and crossing the bound- 
ary between the stable and the unstable regions. By trac- 
ing a single dimensional trajectory between known criti- 
cal points in the infinite dimensional coordinate space we 
can learn something of the higher-dimensional structure. 

FigurcrLH illustrates the topographical variation in the 
energy surface when the unstable region is approached 
from below as an is increased along the 7 = 0.01 line. 
Here, the critical points (including the unstable or delo- 
calized solutions) for each (an, 7) pair were evaluated, 
and linear interpolation was used to generate a continu- 
ous path between stable (S), unstable (U) and delocal- 
ized (D) solutions, along a contour of constant M' |37j ]. 

It is clear from this illustration that, as an ap- 
proaches the unstable region, the soliton solution be- 
comes metastable as it corresponds to a local rather 
than absolute minimum in the energy surface. This fact 
does not affect the existence of stable solutions, however, 
which persist until the local minimum is completely elimi- 
nated as the boundary is crossed into the unstable region. 

We can use this analysis to explain the existence of 
GVA solutions which delocalize under evolution despite 
the existence of an exact stationary solution. In such 
cases, the width of the confining local minimum well (in 
state space) is smaller than the perturbation to the ex- 
act solution induced by enforcing the Gaussian ansatz. 
As one might expect, this only happens near a stability 
boundary, where the wells are small. 



VI. RELATION TO PHYSICAL PARAMETERS 

From the practical point of view, the question of in- 
terest is whether a given set of the physical parameters 
m i, Xi Aw, and the total number of particles N 
can support stable atomic-molecular solitons. In order 
to answer this question we must be able to inter-convert 
between the soliton parameters in terms of the dimen- 
sionless variables and the original physical parameters. 
This procedure is not of a trivial matter, and requires 
a self-consistent solution that is able to map a given set 
of values of (mi, x, ku, Aw, N) into a pair of values 
of (7, an), using the time scale to an d the length scale 
do- Depending on whether or not the pair of values of 
(7, an) is inside the soliton existence domain 7 > and 
< an < I/Ti one can then answer the above question 
and find the soliton parameter values in terms of dimen- 
sional variables, using the results of previous sections. 

We recall that the relationship between the parameters 
of the dimensionless system and the original physical pa- 
rameters are as follows: 

Aw = (7 - 4)/2i , (44) 



kii/x 2 = auto, (45) 



1 / h \ D l 2 

Here, the role of the mass m\ is in setting up the length 
scale do = y/hto/2m\, so that the soliton widths are 
found unambiguously once the corresponding time scale 
to is set up self-consistently. This leaves us, in general, 
with four independent variables (x, Aw, N) in the 
physical parameter space and two independent variables 
(7, an) of the dimensionless system. Therefore, in order 
to be able to map the soliton existence domain in the 
(7, an) plane into the physical parameter space, we have 
to restrict ourselves to cases where - out of four physical 
parameters (x, in, Aw, N) - only two can be varied 
independently, while the other pair must be kept fixed. 

Depending on which pair of the physical parameters is 
chosen to be fixed or varied, one can identify six different 
cases where the solution of the problem can be found un- 
ambiguously. As an example we consider the case where 
the fixed pair of the parameters are the couplings \ an d 
Kn, while the adjustable parameters are the detuning Aw 
and the total number of particles N. This is the most 
physically relevant case, as the detuning and the total 
number of particles are easier to vary experimentally. 

The procedure of mapping the soliton existence do- 
main in the (7, an) plane into the (Aw, AT) parameter 
space consists of solving Eqs. 1)44 (1 -1)46 (1 to firstly iden- 
tify to, using Eq. l|4*5|l . and then finding the respective 
values of Au and N from the remaining two equations. 
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FIG. 14: Variation in N and Aw along the (an, 7) lines shown 
in Fig. Q The values of the couplings x an d "u are held 
constant: x = 10" 6 m 3/2 /s and « = 4.96 x 10~ 17 m 3 /s. 



As an intermediate step, this involves the evaluation of 
N' = N' (7, an) using Eq. 1(151 ). where the soliton pa- 
rameters a, b, A and B are found from the GVA solutions 
for the (7, an) pair in question. Figure 1711 demonstrates 
this mapping for the values of K11 and \ typical of a 87 Rb 
BEC experiments. 

Similar mapping can easily be constructed in other 
cases, where depending on the choice of the fixed pair 
of the physical parameters the sequential order of solv- 
ing Eqs. ^44 (1 -1(46 )1 will vary. In all cases, the initial step 
should consist of identifying the value of the 'dummy' 
parameter to using one of the equations lf4"^l) - (|4*6T) . and 
then eliminating it in favor of the remaining pair of the 
physical parameters in question. 



VII. SUMMARY 

To summarize, we have applied both Gaussian and ex- 
ponential variational approximations to the problem of 
identifying 3D soliton solutions to parametrically cou- 
pled dilute atomic and molecular Bose condensates, with 
atomic self-interaction present. The soliton existence do- 
main has been investigated, and for 7 > found to be 
defined by < an < I/7 in both cases. 

A detailed numerical study of the dynamical behaviour 
of the Gaussian ansatz solutions has revealed that, as a 
rule of thumb, localized propagation of GVA solutions 
indicates the existence of true (and stable) stationary so- 
lutions. Furthermore, in regions where the GVA propa- 
gation is strongly oscillatory, the corresponding exponen- 
tial ansatz solutions have been found to possess a lower 
energy, and propagate almost without oscillations. This 
indicates that systems for which the energy of molecu- 
lar formation is large and negative and which have weak 
atomic self-interaction possess stationary solutions better 
approximated by the exponential rather than Gaussian 
ansatz. Finally, we have identified an anomalous insta- 
bility 'island' in the (an, 7) parameter space, where the 
Gaussian solitons were dynamically unstable and where 
no exact stationary solutions were numerically found. 

The results obtained give the precise conditions under 
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which 3D coupled atomic-molecular BEC solitons can 
form. This is done in terms of the dimensionless parame- 
ters (an, 7) that originate from physical parameters de- 
termined by the atom-molecule coupling, atom-atom re- 
pulsive s -wave scattering, and the energy detuning be- 
tween the atomic and molecular fields. The total num- 
ber of particles N in the system is incorporated in the 
analysis self-consistently, via scaling with respect to the 
time scale to or equivalently with respect to the choice 
of the energy origin. We have shown how one can in 
principle map the physical parameter space into the pa- 
rameter space defined by the dimensionless parameters 
(an, 7), and hence identify whether a particular set of 
physical parameter values, (x, K\\, Aw, N), lays within 
the identified soliton existence and stability domain in 
the (an, 7)-plane. 

Beyond the theoretical interest, the realization of this 
type of coupled BEC solitons could provide a way to 
stabilize the otherwise diverging output of an atom- 
molecular laser, providing a technique for delivering a 
localized, intense and coherent atomic-molecular field to 
a target or a detector. 

The steady-state and variational solutions found here 
are also applicable to generalizations of the Schrodinger- 
Newton equation. 
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FIG. 15: The variation in an as a function of b, for (a) 7 = 
and (b) 7 = 10 s , illustrating that an is monotonically 
decreasing function of b on < b < b q . 
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APPENDIX A 

Standard algebraic analysis of Eqs. (|36[1 - 1)39(1 for pos- 
sible solutions, for a given pair of an > and 7, gives 
the following results. First of all, it is easy to see that 
for 6 > and B > 0, the requirement a > and A 2 > 
(so that A 2 B > too) can be fulfilled only if 7 + b > 0. 
From this, it also follows that 3b + 7 > 0, 56 + 7 > 0, 
and b — 2a < 0. Next, one can easily see that B > if 
a — 1 < 0. This last inequality can be solved in terms of 
6, giving the result that b must satisfy b < b±, where 



1 r 



/,, = - [(2 - 7) + ^(2-7)2+407 



(47) 



corresponds to the positive valued solution of the 
quadratic equation for 6 which follows from a — 1 = 0. 
Thus, as long as < b < b\ (a — 1 < 0) and 7 + b > 0, we 
satisfy the requirements that a > 0, B > 0, and A 2 > 0. 

Analysis of the remaining requirement of an > is 
now reduced to the solution of 4a 2 — ab — b < 0, within 
the region < b < b\ (where a — 1 < 0) and for 7 + 6 > 0. 



Here one has to distinguish and consider separately two 
cases corresponding to 7 < and 7 > 0. 

For 7 < the analysis is quite complicated and requires 
a pure numerical investigation. The overall conclusion is 
that the soliton existence domain is now restricted to a 
very small region which is not of major physical interest. 
For example, for an = 0, this domain is limited to the 
values of 7 within a narrow interval —0.0074 < 7 < 0. 
The size of this interval decreases with increasing an and 
approaches — 1/an < 7 < as an — > 00. In terms of 
the original (physical) phase mismatch parameter Aw, a 
physically interesting and important region corresponds 
to the values of Aw ~ which we note correspond to 
7 ~ 4, while 7 < corresponds to very large and negative 
detuning Alu. For this reason, in the remaining of the 
paper we only treat the case of 7 > 0. 

For 7 > we proceed with the analytic treatment as 
follows. Substituting a from Eq. I|3t)|l . one can rewrite 
the inequality 4a 2 — ab — b < (equivalent to an > 0) 
in the following form 



456 2 + (14 7 - 2)6 + 7(7 - 4) < 



2 7 2 



(48) 



Taking here the equality sign, and considering the left 
and the right hand sides as functions of 6, gives an equa- 
tion that can always be solved graphically, with the result 
that for 7 > there always exists one and only one real 
positive solution for 6, which we define via 6 . Conse- 
quently, the above inequality and therefore an > is 
satisfied if < 6 < 6q. 
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The next step in our analysis consists of showing that 
b < b±, thus restricting the soliton existence region to 
< b < bo from which it follows that for 7 > the 
requirements a > 0, an > 0, B > 0, and A 2 > are 
satisfied simultaneously. The proof of bo < b\ is ac- 
complished by first noting that b — b\ (corresponding to 
a — 1 = 0) is the pole of an, and that an is a continuous 
single- valued function of b within the interval < b < b\ . 
The discontinuity at b = b\ is such that an — > — 00 as 
6 — > 61 — and an — > +00 as + 0. In addition, 

lim^o^n = I/7 > 0, and since bo is the only positive 
root of the cubic equation l|48ll . we conclude that the 
crossover (b — bo) of an from positive to negative values 
can only take place at b < b\ , which implies that 60 < &i • 

Finally, one can show graphically (see Fig. I15J) . that 
within the interval < b < bo and for 7 > 0, an is a 
monotonically decreasing function of b. It approaches its 
maximum an — * 1/7 as b — ► 0, and an = at b = 
b . This implies that there is one-to-one correspondence 
between the soliton existence domain < b < bo and the 
interval < an < 1/7, for 7 > 0. In other words, the 
interval < an < 1/7 (7 > 0) can be equally regarded 
as the soliton existence domain in the parameter space 
("11,7)- 



APPENDIX B 

Here we analyze the set of algebraic equations (1401 )- 
()43|) for possible solutions for a given pair of an > and 
7, with p, q, P and Q all being real and positive. The 
analysis is very similar to the one carried for the GVA 
solutions, and we summarize it as follows. 

For q > (as required), in order that p > it is clearly 
necessary that q 2 + 37 > 0, since in this case we also have 
5q 2 + 37 > 0. In addition, for P 2 to be positive one has to 
have q 2 + 7 > 0, provided that Q > (as required). The 
condition q 2 + 37 > then leads to the requirement that 
p 2 — 3 < for Q to be indeed positive, since q 2 + 3j > 
and q > imply that q — 2p < 0. 

Substituting the expression for p into p 2 — 3, the re- 
quirement p 2 — 3 < can be written in the form 

25<j 6 + (30 7 - 12)q 4 + (9 7 2 - 72~/)q 2 - IO87 2 < 0. (49) 

It is clear that, as the left-hand side of the inequality 
is a polynomial of even order, it can be satisfied if its 
largest real root, q\, is positive. By performing a numer- 
ical analysis of this polynomial, we find that real positive 
solutions exist for all 7, and that q\ — ► 2a/3 as 7 — * 00. As 
in the case of the GVA solutions, we restrict ourselves to 
the physically interesting subspace of 7 > 0. In this case, 
there always exist one and only one real positive solution 
<7i to the above polynomial. Thus, provided < q < q\ 
(i.e. p 2 — 3 < 0) and q 2 + 3j > 0, the requirements that 
p > 0, Q > and P > are met. 



Turning to the remaining equation for an, Eq. I|41|l . 
with an > 0, we next find that we must have 4p 3 — 
q(p 2 + 3) < (in conjunction with p 2 — 3 < 0) on the 
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FIG. 16: The variation in a\\ as a function of q, for (a) 7 = 
and (b) 7 = 10 5 . Note that go < Qi and that qo,qi — > 2y3 
for large 7. 



domain < q < q\ . Here, we have also taken into account 
the fact that the term (q 2 + 7) is always positive, for 
7 > 0. Substituting the expression for p from Eq. 
into 4p 3 — q(p 2 + 3) < we find that this inequality is 
equivalent to 



241q 8 + (3457 - 12)g 6 + 7 (171 7 - 108)q 4 
+7 2 (27 7 - 324)g 2 - 324 7 3 < , 



(50) 



and can be satisfied if the largest real root qo to the 
polynomial in the left-hand side is positive. We again 
employ numerics to find that for 7 > there exist only 
one real and positive root qo, i.e. the above inequality 
is satisfied for q on < q < qo- In the limit of large 7, 
go^2\/3. 

Finally, as for the GVA solutions, we can show nu- 
merically (see Fig. I16[l that an is a monotonically de- 
creasing function of q on < q < qo, and that qo < q\. 
It approaches its maximum value I/7 as q — > 0, and 
an = at q = qo- This implies that there exists one- 
to-one correspondence between the intervals < q < qo 
and < an < 1/7 (7 > 0), so that the later can equally 
be regarded as the soliton existence domain for the EVA 
solutions on the parameter space (an, 7). 
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